
####################################################################################################
#klambda is a function to add dependence into the exponential draws via dependence
#in risk by event. The form of the event dependence is k*lambda_00
###################################################################################################

klambdaexp <- function(risk, k, rho=1){
  n <- length(risk)
  timesbyk<-matrix(0,n,k)
  if(rho==1) {
    timesbyk<-matrix(rexp(n*k,rep(risk,k)),nrow=n)
  }
  else {
    for(i in 1:k){
      timesbyk[,i]<-rexp(n,risk+((i-1)*risk))
    }
  }
  return(timesbyk)
} 

klambdaweibull <- function(risk, shape, k, rho=1){
  n <- length(risk)
  timesbyk<-matrix(0,n,k)
  if(rho==1) {
    timesbyk<-matrix((-log(runif(n*k))/risk)^(1/shape),nrow=n)
  }
  else {
    for(i in 1:k){
      timesbyk[,i]<-(-log(runif(n))/(risk+((i-1)*risk)))^(1/shape)
    }
  }
  return(timesbyk)
}

klambdaweibullelapsed <- function(risk, shape, k, rho=1){
  n <- length(risk)
  timesbyk<-matrix(0,n,k)
  if (rho==1){
    for(i in 1:k){
      if (i==1){
        timesbyk[,i] <- (-log(runif(n))/risk)^(1/shape)
      } else{
        maxunif <- exp(-1*risk*(timesbyk[,i-1]^shape))
        timesbyk[,i] <- (-log(runif(n,0,maxunif))/risk)^(1/shape) 
      }
    }
  }
  else {
    for(i in 1:k){
      if (i==1){
        timesbyk[,i] <- (-log(runif(n))/risk)^(1/shape)
      } else{
        newrisk <- risk+(i-1)*risk
        maxunif <- exp(-1*newrisk*(timesbyk[,i-1]^shape))
        timesbyk[,i] <- (-log(runif(n,0,maxunif))/newrisk)^(1/shape) 
      }
    }
  }
  return(timesbyk)
}
#debug(klambdaweibullelapsed)
#klambdaweibullelapsed(c(rep(exp(0),10),rep(exp(-1),10)),.8,10,1)

klambdanegweibull <- function(risk, shape, k, rho=1){
  n <- length(risk)
  timesbyk<-matrix(0,n,k)
  if(rho==1) {
    timesbyk<-matrix((-log(runif(n*k))/risk)^(1/shape),nrow=n)
  }
  else {
    for(i in 1:k){
      timesbyk[,i]<-(-log(runif(n))/(risk-((i-1)*risk)))^(1/shape)
    }
  }
  return(timesbyk)
} 

klambdagompertz <- function(risk, shape, k, rho=1){
  n <- length(risk)
  timesbyk<-matrix(0,n,k)
  if(rho==1){
    timesbyk<-matrix((1/shape)*log(1-(shape*log(runif(n*k))/risk)),nrow=n)
  }
  else {
    for(i in 1:k){
      timesbyk[,i]<-matrix((1/shape)*log(1-(shape*log(runif(n))/(risk+(i-1)*risk))),nrow=n)
    }
  }
  return(timesbyk)
}


covarunif <- function(len,min=-2,maxi=2) {
  r2unif<-runif(len,min,maxi)
  return(r2unif)
}

covarnorm <- function(len,mean=0,std=1) {
  none<-vector("integer",len)
  if(std==0){
    renorm<-none
  }
  else { 
    renorm<-rnorm(len,mean,std)
  }
  return(renorm)
}


covargamma <- function(len,theta=1) {
  none<-vector("integer",len)
  if(theta==0){
    regamma<-1
  }
  else { 
    regamma<-rgamma(len,shape=(1/theta),scale=theta)
  }
  return(regamma)
}


simdataexp <- function(nn, maxev, maxt, p, rxe, base, covar, times, rho=1, min=0, maxi=1) {
  n <- floor(nn/2)
  xx <- covar(nn, min, maxi)
  rx <- c(rep(1, n), rep(0, n))
  risk <- xx
  risk[1:n] <- risk[1:n] + rep(rxe, round(p * n))
  rxgrp <- c(rep(1:length(p), round(p * n)), rep(0, n))
  gapstart <- rep(0, nn)
  risk <- exp(risk) * base 
  temp <- times(risk, maxev, rho)
  temp2 <- apply(temp, 1, cumsum)
  temp3 <- rbind(0, temp2[-maxev, ])
  keep <- (temp3 < maxt)
  stat <- 1 * (temp2 < maxt)
  temp2 <- ifelse(temp2 < maxt, temp2, maxt)
  start2 <- c(temp3[keep])
  stop2 <- c(temp2[keep])
  gapstop <- stop2 - start2
  enum0 <- (rep(1:maxev - 1, nn))[keep]
  strata <- ifelse(enum0 == 0, 0, ifelse(stop2 < maxt, enum0, 
                     ifelse(stop2 == maxt & enum0 == 1, enum0, enum0 - 1)))
  td <- data.frame(id = (rep(1:nn, rep(maxev, nn)))[keep], 
                   start = c(temp3[keep]), stop = c(temp2[keep]),
                   status = stat[keep],
                   rx = (rep(rx, rep(maxev, nn)))[keep],
                   x = (rep(xx, rep(maxev,nn)))[keep],
                   enum = (rep(1:maxev, nn))[keep],
                   rxgrp = (rep(rxgrp, rep(maxev, nn)))[keep],
                   gapstart = (rep(gapstart,rep(maxev, nn)))[keep],
                   gapstop = c(gapstop),
                   strata = c(strata))
  for (i in 1:maxev) td[[paste("rx", i, sep = "")]] <- td$rx * 
    (td$enum == i)
  if (maxev > 2) {
    for (i in 2:(maxev - 1)) td[[paste("rx", i, "p", sep = "")]] <- td$rx * 
      (td$enum >= i)
  }
  td
}


simdatashape <- function(nn, shape, maxev, maxt, p, rxe, base, covar, times, rho=1, min=0, maxi=1) {
  n <- floor(nn/2)
  xx <- covar(nn, min, maxi)
  rx <- c(rep(1, n), rep(0, n))
  risk <- xx
  risk[1:n] <- risk[1:n] + rep(rxe, round(p * n))
  rxgrp <- c(rep(1:length(p), round(p * n)), rep(0, n))
  gapstart <- rep(0, nn)
  risk <- exp(risk) * base
  temp <- times(risk, shape, maxev, rho)
  temp2 <- apply(temp, 1, cumsum)
  temp3 <- rbind(0, temp2[-maxev, ])
  keep <- (temp3 < maxt)
  stat <- 1 * (temp2 < maxt)
  temp2 <- ifelse(temp2 < maxt, temp2, maxt)
  start2 <- c(temp3[keep])
  stop2 <- c(temp2[keep])
  gapstop <- stop2 - start2
  enum0 <- (rep(1:maxev - 1, nn))[keep]
  strata <- ifelse(enum0 == 0, 0, ifelse(stop2 < maxt, enum0, 
                     ifelse(stop2 == maxt & enum0 == 1, enum0, enum0 - 1)))
  td <- data.frame(id = (rep(1:nn, rep(maxev, nn)))[keep], 
                   start = c(temp3[keep]), stop = c(temp2[keep]),
                   status = stat[keep],
                   rx = (rep(rx, rep(maxev, nn)))[keep],
                   x = (rep(xx, rep(maxev,nn)))[keep],
                   enum = (rep(1:maxev, nn))[keep],
                   rxgrp = (rep(rxgrp, rep(maxev, nn)))[keep],
                   gapstart = (rep(gapstart,rep(maxev, nn)))[keep],
                   gapstop = c(gapstop),
                   strata = c(strata))
  for (i in 1:maxev) td[[paste("rx", i, sep = "")]] <- td$rx * 
    (td$enum == i)
  if (maxev > 2) {
    for (i in 2:(maxev - 1)) td[[paste("rx", i, "p", sep = "")]] <- td$rx * 
      (td$enum >= i)
  }
  td
}


simdatashapegamma <- function(nn, shape, maxev, maxt, p, rxe, base, times, rho=1, theta=0) {
  n <- floor(nn/2)
  xx <- covarnorm(nn, mean=0,std=0)  #for this model, this does nothing but create a vector of 0's length n
  rx <- c(rep(1, n), rep(0, n))
  risk <- xx
  risk[1:n] <- risk[1:n] + rep(rxe, round(p * n))
  rxgrp <- c(rep(1:length(p), round(p * n)), rep(0, n))
  gapstart <- rep(0, nn)
  frail <- covargamma(nn, theta=theta)
  risk <- exp(risk) * base * frail
  temp <- times(risk, shape, maxev, rho)
  temp2 <- apply(temp, 1, cumsum)
  temp3 <- rbind(0, temp2[-maxev, ])
  keep <- (temp3 < maxt)
  stat <- 1 * (temp2 < maxt)
  temp2 <- ifelse(temp2 < maxt, temp2, maxt)
  start2 <- c(temp3[keep])
  stop2 <- c(temp2[keep])
  gapstop <- stop2 - start2
  enum0 <- (rep(1:maxev - 1, nn))[keep]
  strata <- ifelse(enum0 == 0, 0, ifelse(stop2 < maxt, enum0, 
                     ifelse(stop2 == maxt & enum0 == 1, enum0, enum0 - 1)))
  td <- data.frame(id = (rep(1:nn, rep(maxev, nn)))[keep], 
                   start = c(temp3[keep]), stop = c(temp2[keep]),
                   status = stat[keep],
                   rx = (rep(rx, rep(maxev, nn)))[keep],
                   x = (rep(xx, rep(maxev,nn)))[keep],
                   enum = (rep(1:maxev, nn))[keep],
                   rxgrp = (rep(rxgrp, rep(maxev, nn)))[keep],
                   gapstart = (rep(gapstart,rep(maxev, nn)))[keep],
                   gapstop = c(gapstop),
                   strata = c(strata))
  for (i in 1:maxev) td[[paste("rx", i, sep = "")]] <- td$rx * 
    (td$enum == i)
  if (maxev > 2) {
    for (i in 2:(maxev - 1)) td[[paste("rx", i, "p", sep = "")]] <- td$rx * 
      (td$enum >= i)
  }
  td
}



simdatashapegamma.random <- function(nn, shape, maxev, censorrate, maxt, p, rxe, base, times, rho=1, theta=0) {
  n <- floor(nn/2)
  xx <- covarnorm(nn, mean=0,std=0)  #for this model, this does nothing but create a vector of 0's length n
  rx <- c(rep(1, n), rep(0, n))
  risk <- xx
  risk[1:n] <- risk[1:n] + rep(rxe, round(p * n))
  rxgrp <- c(rep(1:length(p), round(p * n)), rep(0, n))
  gapstart <- rep(0, nn)
  frail <- covargamma(nn, theta=theta)
  risk <- exp(risk) * base * frail
  temp <- times(risk, shape, maxev, rho) # creates matrix of random durations (nn x maxev)
  temp2 <- apply(temp, 1, cumsum) # creates stop-time matrix (maxev x nn)
  temp3 <- rbind(0, temp2[-maxev, ]) # creates start-time maxtrix (maxev x nn)

  keep <- (temp3 < maxt)
  stat <- 1 * (temp2 < maxt)
  temp2 <- ifelse(temp2 < maxt, temp2, maxt)
  
  if (censorrate==0){    
    ## No Censoring
    start2 <- c(temp3[keep])
    stop2 <- c(temp2[keep])
    gapstop <- stop2 - start2
    enum0 <- (rep(1:maxev - 1, nn))[keep]
    strata <- ifelse(enum0 == 0, 0, ifelse(stat[keep]==1,enum0, ifelse(stat[keep]==0 & enum0 == 1, enum0, enum0 - 1)))
  } else {
    ## Random Censoring
    cases <- 1:nn  # index all cases, may need to fine-tune censorrate precision
    while (sum(stat[keep])/sum(1*keep) >= (1-censorrate)){
      case.omit <- sample(cases,1)  # Allow for resampling
      event.omit <- ceiling(runif(1,min=0,max=maxev)) # Select event (out of max) to censor
      if (stat[event.omit,case.omit]==1){  # If observation isn't already censored, make it so
        stat[event.omit:maxev,case.omit] <- 0  # Indicate censoring
        temp2[event.omit:maxev,case.omit] <- runif(1,min=temp3[event.omit,case.omit],max=temp2[event.omit,case.omit])  # Change stop time to random censor point
        if (event.omit<10){
          keep[(event.omit+1):maxev,case.omit] <- as.logical(0) # Remove events after censored observation
        }
      }
    }
    start2 <- c(temp3[keep])
    stop2 <- c(temp2[keep])
    gapstop <- stop2 - start2
    enum0 <- (rep(1:maxev - 1, nn))[keep]
    strata <- ifelse(enum0 == 0, 0, ifelse(stat[keep]==1,enum0, ifelse(stat[keep]==0 & enum0 == 1, enum0, enum0 - 1)))
  }

  td <- data.frame(id = (rep(1:nn, rep(maxev, nn)))[keep], 
                   start = c(temp3[keep]), stop = c(temp2[keep]),
                   status = stat[keep],
                   rx = (rep(rx, rep(maxev, nn)))[keep],
                   x = (rep(xx, rep(maxev,nn)))[keep],
                   enum = (rep(1:maxev, nn))[keep],
                   rxgrp = (rep(rxgrp, rep(maxev, nn)))[keep],
                   gapstart = (rep(gapstart,rep(maxev, nn)))[keep],
                   gapstop = c(gapstop),
                   strata = c(strata))
  for (i in 1:maxev) td[[paste("rx", i, sep = "")]] <- td$rx * 
    (td$enum == i)
  if (maxev > 2) {
    for (i in 2:(maxev - 1)) td[[paste("rx", i, "p", sep = "")]] <- td$rx * 
      (td$enum >= i)
  }
  td
}

elapsedsimdatashapegamma <- function(nn, shape, maxev, maxt, p, rxe, base, times, rho=1, theta=0) {
  n <- floor(nn/2)
  xx <- covarnorm(nn, mean=0,std=0)  #for this model, this does nothing but create a vector of 0's length n
  rx <- c(rep(1, n), rep(0, n))
  risk <- xx
  risk[1:n] <- risk[1:n] + rep(rxe, round(p * n))
  rxgrp <- c(rep(1:length(p), round(p * n)), rep(0, n))
  gapstart <- rep(0, nn)
  frail <- covargamma(nn, theta=theta)
  risk <- exp(risk) * base * frail
  temp2 <- t(times(risk, shape, maxev, rho))
  temp3 <- rbind(0, temp2[-maxev, ])
  keep <- (temp3 < maxt)
  stat <- 1 * (temp2 < maxt)
  temp2 <- ifelse(temp2 < maxt, temp2, maxt)
  start2 <- c(temp3[keep])
  stop2 <- c(temp2[keep])
  gapstop <- stop2 - start2
  enum0 <- (rep(1:maxev - 1, nn))[keep]
  strata <- ifelse(enum0 == 0, 0, ifelse(stop2 < maxt, enum0, 
                     ifelse(stop2 == maxt & enum0 == 1, enum0, enum0 - 1)))
  td <- data.frame(id = (rep(1:nn, rep(maxev, nn)))[keep], 
                   start = c(temp3[keep]), stop = c(temp2[keep]),
                   status = stat[keep],
                   rx = (rep(rx, rep(maxev, nn)))[keep],
                   x = (rep(xx, rep(maxev,nn)))[keep],
                   enum = (rep(1:maxev, nn))[keep],
                   rxgrp = (rep(rxgrp, rep(maxev, nn)))[keep],
                   gapstart = (rep(gapstart,rep(maxev, nn)))[keep],
                   gapstop = c(gapstop),
                   strata = c(strata))
  for (i in 1:maxev) td[[paste("rx", i, sep = "")]] <- td$rx * 
    (td$enum == i)
  if (maxev > 2) {
    for (i in 2:(maxev - 1)) td[[paste("rx", i, "p", sep = "")]] <- td$rx * 
      (td$enum >= i)
  }
  td
}

